Quantum simulation of manybody effects in steady-state nonequilibrium: 
electron-phonon coupled quantum dots 
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We develop a mapping of quantum steady-state nonequilibrium to an effective equilibrium and 
solve the problem using a quantum simulation technique. A systematic implementation of the 
nonequilibrium boundary condition in steady-state is made in the electronic transport on quantum 
dot structures. This formulation of quantum manybody problem in nonequilibrium enables the use of 
existing numerical quantum manybody techniques. The algorithm coherently demonstrates various 
transport behaviors from phonon-dephasing to I — V staircase and phonon- assisted tunneling. 

PACS numbers: 73.63.Kv, 72.10.Bg, 72.10.Di 
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Electronic transport in nanoscale devices has recently 
attracted much attention for its tremendous impact in 
modern electronics. In nanoscale devices, fundamental 
understanding of quantum physics holds the key to a 
breakthrough in the field. One of the central issues is to 
understand the increasing interplay of quantum many- 
body effect and the high-bias nonequilibrium. 

Despite the intense research in the nanoscale transport, 
most of the theoretical tools are based on the Green func- 
tion technique 0, 0, 0, , and powerful numerical equi- 
librium techniques, which have been essential in quan- 
tum manybody theories, have not been useful due to 
the lack of nonequilibrium formulation. Application of 
such well-established numerical techniques, such as quan- 
tum Monte Carlo (QMC) method 0, to nonequilibrium 
with full quantum manybody effects will provide a break- 
through in the field. The goal of this work is to formulate 
a general scheme of combining the nonequilibrium the- 
ory with numerical quantum manybody techniques via a 
mapping of a nonequilibrium system to an effective equi- 
librium. This method applied in electron-phonon coupled 
quantum dot systems demonstrates various transport be- 
havior such as the phonon dephasing, phonon-assisted 
tunneling, I — V staircase in a single framework. The 
formulation shown is applicable to other numerical tech- 
niques, such as nemerical renormalization group, exact 
diagonalization, etc. 

Conventional diagrammatic theory is founded in adia- 
batic turn-on of interaction. Here we begin from a very 
different point that a steady-state nonequilibrium pos- 
sesses the time-invariance, where the transient states are 
damped out. With this observation, Hershfield [f| has 
shown the existance of an effective Hamiltonian which 
builds the nonequilibrium ensemble in steady-state. The 
ensemble can be represented in terms of a modified par- 
tition function Z„ 



-'noncq for nonequilibrium as 
= Tre -P(H-Y) with p-i/T, 
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FIG. 1: (a) Schematic diagram of a quantum transport de- 
vice. The source and drain reservoirs are biased by a chemical 
potential difference 3>. The quantum dot (QD) level is set in 
the middle for simplicity, (b) Each reservoir is mapped to 
a single site and a new reservoir. The three discrete sites 
for the QD and L- and .R-sites are treated in fully quantum 
mechanical manner. 



reservoirs and quantum dot (QD) (see FIG. [I] (a)). The 
additional operator Y accounts for the boundary condi- 
tion of chemical potential difference in the two reservoirs, 
given by the bias <&. Once the bias operator Y is con- 
structed, one can apply the usual equilibrium QMC tech- 
nique [a] to -f/noneq = H — Y as if -ff n oneq describes an 
equilibrium. To obtain I—V characteristics, for instance, 
one simply needs to evaluate the expectation value of the 
current operator I, 



Tr Ie-I 3 ^-^ 
Tr e-P(H-Y) 
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where H is the usual Hamiltonian operator for the whole 
system including the source (left, L), drain (right, R) 



(2) 

where c^ Lka {chka) is the electron creation (annihilation) 
operator of spin <r, the fc-th continuum state in the L- 
reservoir, d\ (da-) the operator for QD state, and th the 
hopping integral (with the volume factor f2) between the 
QD and the i-reservoir (see FIG.^a)). 

The central challenge is to construct the bias opera- 
tor Y . Once a nonequilibrium ensemble is established 
via e~P^ H ~ Y \ the subsequent motion is governed by the 
time evolution operator e~ tHt . Therefore the bias opera- 
tor Y must satisfy the commutation relation [H, Y] = 0. 
Hershfield 61 has shown that Y can be written in terms 



2 



of the scattering state operator ip^ akcr (a = L,R) as Q 



Y = 



2~^ 



(3) 



The scattering state creation operator i$ aklJ satisfies the 
Lippman-Schwinger equation in the operator form, 

[H, V'LJ = Confer + *»?(^fca - c akJ> ( 4 ) 

with e Q fe the asymptotic continuum energy and r\ an in- 
finitesimal convergence factor. 

There are two main issues to implement the above 
nonequilibrium ensemble. First, to systematically con- 
struct Y, the scattering state is expanded in terms of 
interaction as ^ aka = £ n V>L CT) „ H satisfying 

\V^lka,n-l] = ( e «fc +*V)^aka, n ~ V&to.J. ( 5 ) 

where the total Hamiltonian H = Hq + V. Hq is the 
non-interacting and V the interacting part of the Hamil- 
tonian. The second and more serious issue is the non- 
locality of the operator Y. Since the one-particle eigen- 
state of the total Hamiltonian is delocalizcd in space, 
the bi-products ip^^aka in Y produce non-local terms 
even if the model Hamiltonian H is given with local in- 
teractions. Despite the great potential of Hershfield's pa- 
per , only limited developments @, Q have been made 
so far due to such difficulties. 

To make the following discussions concrete, let us con- 
sider an example: electron-phonon (el-ph) coupled QD 
connected to two L, i?-reservoirs. Spinless Hamiltonian 
for the system reads H = H + V = Hf + H% h + V, 

Hq = ^ £qfcc^ fc c a fc + € d Sd + (d^Cgk + h.c.) 



ak ak 

Hf = \{P 2 + ^W) , V = ap(dU - (dU)), (6) 



where ip is phonon amplitude, p its conjugate momen- 
tum, w p h phonon frequency, and a(= g^J2uj p h) the el-ph 
coupling constant. As shown in FIG. [Ja), we model 
the continuum states to be shifted with the bias, i.e., 
(Lk = ?k + $/2 and e Rk = e k — $/2. We set the particle- 
hole symmetric QD level (td = 0) for simplicity. 

Here we propose that we capture important many- 
body effects by isolating the leading order non-local in- 
teractions as follows. We make two observations in 
H — Y . First, the QD level is modulated by phonons 
as ed(<p) = £d + otip. The fluctuation of the QD level 
results in dephasing of the current. Second, the hop- 
ping out of the QD to the continuum is effectively mod- 
ified as e d (p) is driven in and out of resonance with re- 
spect to the reservoir chemical potentials. As will be 
discussed later, fluctuations in the hopping integral re- 
sult in phonon satellites in current. With these observa- 
tions, we isolate the first path of hopping from the QD 



by rewriting the each continuum as a fictitious site and a 
new continuum, as depicted in FIG.^b). We then treat 
the coupling between the three sites fully quantum me- 
chanically and make mean-field approximations to terms 
involving the continuum states. This procedure can be 
systematically improved by inserting additional fictitious 
reservoir sites to a better non-local approximation. 

We have used a semi-circular density of state (DOS) for 
the continua, N(e a ) = — e^/irD 2 , with the Fermi 

energy ej? = D. The semi-circular DOS is particularly 
useful since, with the intra-reservoir hopping to = D/2 
(see FIG. [Q, the resulting new DOS is identical to the 
original DOS. In the non-interacting limit (V = 0), the 
scattering state can be readily obtained by expanding 
ijy a1t in terms of the basis states in Eq. Q as 



to 



-Lk 



^29°L,p{tLk)c} r 



(7) 



where /3 = Lk' , L,d,R, Rk' with L, R denoting the ficti- 
tious reservoir sites, and Lk' , Rk' new continuum states. 

is the same as d^. g°^(e) are retarded Green func- 
tions propagating from state L to j3. ipjfa are similarly 
obtained. From Eq. J2J), one can recover the Landauer- 
Biittiker formula in the non- interacting limit. 

We expand the scattering states up to harmonic 
el-ph coupling via Eq. JHJ. With [V'jV'iko] = 
(t /VU)aipg° Ld (e Lk )d' i and ip' ak 1 expanded as 



,LkJ 



LkJ 



(8) 



one obtains explicit expressions for ap, bp after a straight- 
forward calculation. Up to the linear order of (p,p) we 
obtain Y = % + Y with Y = $/2 E fe (V>LoV>£fco - 
^ma^Rko), Yi = $/2£ fc (V>L V>z,fci + ^'LiV'ifcQ - 
^Rko^ 'R kl ~ V'ljfciV'-Rfco)- Now we make a mean-field ap- 
proximations (ip — 0) on the terms which involve the 
continuum states Lk', Rk'. Projected onto the discrete 
basis (L,d,R), Y\ reads 

Yi = ipA + pB = & A Pi + P B Pi) 4 C 7> ( 9 ) 

0,j=L,d,R 

where the 3x3 matrix A is written as 

A M = ^J2[4 k (9°L^Lk)r+g%(e L k)(a^y 

k 

-af {g%(e Rk ))*-g° R ^e Rk ) (of )*] , (10) 

and B is given by replacing a^ k by bp k . 

We now sample the ensemble by treating H noncq as in 
conventional QMC We express the Boltzmann fac- 
tor into the Trotter breakup [5|, Tr exp(— /3iJ nonoq ) = 



Tr[exp(-Ari? nonoq )] A ' with (3 = NAt. The tl 
crete sites (L,d,R) are considered as a general 
purity and the Hirsch-Fye algorithm 01 is cm F 
integrate out the continuum states which are 
rated in the non-interacting Green function 
the Matsubara imaginary frequency uj n — (2n + '. 

no c \ _ 9%{^k) [s° 7 (e Q fc)]* ^ < 

G M ^ n )- ^ +( _ )a$/2 +2.- 

a=L,R;k V / m 

where the possible localized states (with the m-t 
E m and amplitudes c™ for the (3 state) are ta 
account. (-) L = 1, (—) R = -1. The localize 
become important in the narrow bandwidth lim 
The main difference from the usual equilibriu 
is the p-term in Eq. 0. Integrating out the c 
momentum p inside time-slices of the Trotter- 
leads to [HI 

(12) 

where B T is the electronic operator at the imaginary time 
r. Due to the second term B T {dip/dr), a temporally lo- 
cal update of the phonon field tp T results in changes in 
B t -At, B t , B t +at- The apparent broken time- reversal 
symmetry in the first-order time differentiation (dtp/dr) 
originates from the right-to-left flow of electrons induced 
by the bias $ 12]. The (dtp/dr) term also shows bias- 
induced phonon fluctuations which become crucial in 
phonon-assisted tunneling 0]. At = 1 in the follow- 
ing results. 

I — V curves in the broad band limit (D ^> to pn ) are 
presented in FIG. [3 The unit of energy is chosen so 
that the bare phonon level (oj p h = 0.5) aligns with the 
Fermi energy of the L-reservoir at unit bias ($ = 1). 
The current monotonically increased with $ at all el-ph 
coupling constant g. The current converges toward the 
well-known large-band, large-bias limit (dashed line for 
9 = 0) 

1=^ /^* ' with r a = 7r£iV c (0). (13) 
n i l + J- r 

At zero bias, the QD spectral function Ad(u>) with real 
frequenc y ui (see inset) by using the maximum entropy 
method [15j shows that the phonon satellite peaks are 
clearly separated from the main peak [J^. From this, 
one might have expected that there would be a double- 
step I — V characteristics, while the computed curves do 
not show such distinct features. This suggests that the 
actual decay rate of the phonon satellite under finite bias 
($ w 2wph) is much larger than the zero bias case because 
the phonon level is now close to a resonance with the 
reservoir Fermi energy while at zero bias finite particle- 
hole excitation energies are required for phonons to decay. 
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FIG. 2: Simulated I — V characteristics in the broad band 
limit, D ^> ojph- Parameters are tt, = ta = 0.1, w p h = 

0. 5, D = 8.0 and T = 1/32 = 0.03125. As the electron- 
phonon coupling constant g is increased, the current de- 
creased. Non-interacting phonon level aligns to the reservoir 
Fermi energy at bias $ = 2u v h = 1. The phonon dephasing 
effect is so strong that the phonon satellite peak is invisible at 
all coupling constant while the phonon satellites at zero bias 
(inset) are clearly separated from the main peak. 

However, upon close inspection of the curves, one finds 
that the line shapes change slightly. To guide the eye, a 
tangential line is drawn for g = 0.8 from the zero bias 
where the renormalized coherent current dominates. As 
$ increases, it quickly gives way to the transport which 
involves incoherent phonon excitations. In the broad 
band limit, the dominant role of phonons is the dephasing 
effect from the fluctuation of the QD level td{f) = td+o^P 
via the Yq term, instead of phonon excitations serving as 
well-defined discrete levels. Here the current is rather 
insensitive to the choice of Y\ . 

In the narrow band limit, one obtains a negative differ- 
ential resistance (NDR) behavior ^3]. FIG. shows the 
results when the bandwidth D is gradually decreased. 
Parameters are uj p h = 0.5, = = 0.1 (except in 
(c) as shown) and T — 1/12 = 0.0833. As D becomes 
comparable to <&, a gradual NDR shows up due to the de- 
creasing overlap of L — R DOS. As D is further decreased 
a distinct discrete three-peak feature emerges as shown 
in FIG. |3Jb-c). The peak at near the zero bias is the 
coherent current corresponding to the ballistic transport 
renormalized by the interaction. As D is reduced well be- 
low T, the maximum height of the ballistic current gets 
reduced due to the thermal dephasing. 

The phonon satellite peak at higher bias ($ = 2u pn = 

1, marked as A in the plot) becomes evident in FIG.|3jb). 
The electronic transport in this regime is sequential [lj, 
i.e., electron from the L-reservoir tunnels to the QD by 
emitting a phonon and loses its phase information before 
it tunnels out to the i?-reservoir after emitting another 
phonon. The width of the structure is twice as large as 
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FIG. 3: I — V characteristics for narrow bandwidth limit. 
th — tn — 0.1 (except for (c)), ui ph = 1, g = 0.2 and 
T = 1/12 = 0.0833. (a) The saturated current in the large 
bandwidth limit gradually turns into negative differential re- 
sistance regime as $ becomes comparable to D. (b) Three- 
peak structure emerges as D gets further reduced. The peaks 
correspond to ballistic coherent transport, phonon-assisted 
tunneling (B, $ = uj v h) and sequential off-resonant trans- 
port (A, $ = 2u) p h)- (c) Reduction of tunneling parameter 
th,tn dramatically suppresses the phonon-assisted tunneling 
compared to the sequential current. 
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FIG. 4: (a) Initial and (b) final states in the phonon assisted 
tunneling process. The particle-hole and phonon excitations 
are resonant at the bias $ = uo p h- 



the full bandwidth, 4D, as expected. It is interesting 
that the phonon structures for D = 0.20 have the famil- 
iar stair-case I — V curve, which shows that the phonon 
satellite states act like discrete current channels due to 
the reduced decay paths to continuum, in contrast to the 
large bandwidth limit where the phonon acts more as a 
dephaser than as a discrete level. We emphasize that the 
phonon satellite features A and B are due to the first 
order correction to Y, particularly the p-term in Y\. 

An intriguing feature in the I—V curves is the phonon- 
assisted tunneling marked as B in FIG. As depicted 
in FIG. 01 the tunneling happens by energy-exchange 



between a particle-hole pair across the device and a 
phonon on the QD at the bias <& = ui p h- The curves at 
D = 0.4, 0.2 resemble the experimental results 13]. This 
is a direct consequence of the coherence built in the Y 
operator, Eq. 10, and should be distinguished from the 
sequential tunneling. To demonstrate the point, a calcu- 
lation with reduced tunneling parameters tr, = tji = 0.05 
is shown as a thin line in FIG. Etc). Since the phonon- 
assisted tunneling is through a coherent electron-hole 
state, the tunneling amplitude is proportional to its cross- 
lead overlap thtu, while the amplitude of the sequential 
tunneling is proportional to tx, or tn separately due to 
the uncorrelated sequential tunneling events. The much 
reduced peak B supports the argument. 

We have formulated a quantum simulation algorithm 
for steady-state nonequilibrium and have shown that the 
method reproduced the most important physics in the 
electron-phonon coupled quantum dot system. In com- 
bination with other numerical manybody techniques, the 
formulation of nonequilibrium shown here will provide a 
critical step toward a more complete theory of steady- 
state nonequilibrium. 
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